Agriculture and plant breeding

Images can be used in plant phenotyping to draw inference about many traits:

• Geometric traits (i.e. plant height, leaf area index, lodging, crop canopy cover)
• Canopy spectral texture (spectral features)
• Physiological traits (i.e., chlorophyll, biomass, pigment content, photosynthesis)
• Abiotic/biotic stress indicators (i.e., stomatal conductance, canopy temperature difference, leaf water potential, senescence index)
• Nutrients (nitrogen concentration, protein content)
• Yield

FIELDimageR

FIELDimageR is a R package to analyze images and plant phenotyping, and allows to:

• Crop the image
• Remove soil effect
• Build vegetation indices
• Rotate the image
• Build the plot shapefile
• Extract information for each plot
• Evaluate stand count, canopy percentage, and plant height

Required packages

• FIELDimageR

## Installing 
# install.packages("sp")
# install.packages("raster")
# install.packages("rgdal")
# install.packages("ggplot2")
# install.packages("agricolae")
# install.packages("reshape2")
# install.packages("devtools")
# install.packages("ggrepel")
# install.packages("lme4")
# install.packages("plyr")
# install.packages("DescTools")
# devtools::install_github("filipematias23/FIELDimageR")

## Necessary packages 
library(FIELDimageR)
library(raster)
library(agricolae)
library(reshape2)
library(ggplot2)
library(lme4)
library(plyr)
library(ggrepel)
library(DescTools)

Data (R/FIELDimageR)

• Potato breeding: 140 plots (3x15 feet)
• Population: 97 genotypes
• Six UAV flights: 30 (soil base), 40, 70, and 100 days after planting (DAP)

• Download EX1: RGB orthomosaics.tif MOSAIC and dsm.tif DSM

# EX1 data files: 
MOSAIC<-list.files("./MOSAIC/")
DSM<-list.files("./DSM/")

Basic steps

# Uploading an example mosaic
Test <- stack(paste("./MOSAIC/",MOSAIC[3],sep = ""))
plotRGB(Test)

# Cropping and reducing size 
Test.Crop <- fieldCrop(mosaic = Test)

# Rotating image with theta 2.3:
Test.Rotate<-fieldRotate(Test.Crop, theta = 2.3)
#Test.Rotate<-fieldRotate(Test.Crop, clockwise = F)

# Removing soil (example using HUE):
Test.RemSoil<-fieldMask(Test.Rotate) 

Building the plot shapefile

• Dataset - Field notes - phenotype (Download: ‘EX1_Data.csv’)

# Reading FieldData.csv
Data<-read.csv(file = "EX1_Data.csv",header = T)
Data

# Making the field Map
Map<-fieldMap(fieldPlot = Data$Plot,
              fieldColumn = Data$Column,
              fieldRow = Data$Row,
              decreasing = T)
Map 

• EX1 field design (MAP) should start from bottom to top and left to the right. In other words, plot 1 should be where the blue arrow is showing in the image below.

# Rotating the MAP to correct plots position:
rotate <- function(x) t(apply(x, 2, rev))
Map<-rotate(rotate(Map))
Map

# Building the plot shapefile (ncols = 14 and nrows = 10)
# x11()
plotShape<-fieldShape(mosaic = Test.RemSoil$newMosaic, 
                      ncols = 14, 
                      nrows = 10, 
                      fieldData = Data, 
                      ID = "Plot", 
                      fieldMap = Map)

# Building indices ("NGRDI","GLI")
Test.Indices<- fieldIndex(mosaic = Test.RemSoil$newMosaic, 
                         Red = 1, Green = 2, Blue = 3, 
                         index = c("NGRDI","GLI"), 
                         myIndex = c("(Red-Blue)/Green"))

plot(Test.Indices$myIndex.2)
plot(plotShape$fieldShape,add=T)

Extracting Data

Test.Info<- fieldInfo(mosaic = Test.Indices[[c("NGRDI","GLI","myIndex")]],
                      fieldShape = plotShape$fieldShape,
                      buffer = -0.10,
                      n.core = 3)
Test.Info$fieldShape@data

Estimation of plant height using digital surface model (DSM).

# Uploading files from soil base and vegetative growth:
DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
DSM1 <- stack(paste("./DSM/",DSM[3],sep = ""))
plot(DSM1)

# Cropping the image using the previous shape:
DSM0.Crop <- fieldCrop(mosaic = DSM0,fieldShape = Test.Crop)
DSM1.Crop <- fieldCrop(mosaic = DSM1,fieldShape = Test.Crop)

# Canopy Height Model (CHM):
DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
CHM <- DSM1.Crop-DSM0.R
plot(CHM)

# Rotating the image using the same theta from RGB mosaic:
CHM.Rotate<-fieldRotate(CHM, theta = 2.3)

# Removing the soil using the mask from RGB mosaic:
CHM.RemSoil <- fieldMask(CHM.Rotate, mask = Test.RemSoil$mask)

# Observing EPH profile for 2 ranges:
CHM.Draw <- fieldDraw(mosaic = CHM.Rotate,
                            ndraw = 2)

dev.off()
par(mfrow=c(1,3))

# Profile plot:
plot(x = CHM.Draw$Draw1$drawData$x-100230, 
     y = CHM.Draw$Draw1$drawData$layer, 
     type="l", col="red",lwd=1,ylim=c(0,1),
     xlab="Distance (m)", ylab="EPH (m)")
lines(x = CHM.Draw$Draw2$drawData$x-100230, 
      y = CHM.Draw$Draw2$drawData$layer,
      type="l", col="blue",lwd=1,add=T)

# CHM plot:
plot(CHM.Rotate, col = grey(1:100/100), axes = FALSE, box = FALSE,legend=F)
lines(CHM.Draw$Draw1$drawData$x,CHM.Draw$Draw1$drawData$y, type="l", col="red",lwd=2)
lines(CHM.Draw$Draw2$drawData$x,CHM.Draw$Draw2$drawData$y, type="l", col="blue",lwd=2)

#RGB plot:
plotRGB(Test.Rotate)
lines(CHM.Draw$Draw1$drawData$x,CHM.Draw$Draw1$drawData$y, type="l", col="red",lwd=2)
lines(CHM.Draw$Draw2$drawData$x,CHM.Draw$Draw2$drawData$y, type="l", col="blue",lwd=2)

# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(CHM.RemSoil$newMosaic, 
                 fieldShape = Test.Info$fieldShape, 
                 fun = "mean") #fun="quantile"
EPH$plotValue

Evaluating all mosaics in a loop

DataTotal<-NULL
for(i in 2:length(MOSAIC)){
  EX1 <- stack(paste("./MOSAIC/",MOSAIC[i],sep = ""))
  EX1.Crop <- fieldCrop(mosaic = EX1,
                    fieldShape = Test.Crop,
                    plot = F)
  EX1.Rotate<-fieldRotate(EX1.Crop,
                    theta = 2.3,
                    plot = F)
  EX1.RemSoil<-fieldMask(EX1.Rotate, plot = F)
  EX1.Indices<- fieldIndex(mosaic = EX1.RemSoil$newMosaic, 
                            Red = 1, Green = 2, Blue = 3, 
                            index = c("NGRDI","BGI", "GLI","VARI"), 
                            myIndex = c("(Red-Blue)/Green"),
                            plot = F)
  EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI", "GLI","VARI","myIndex")]],
                        fieldShape = plotShape$fieldShape, 
                        n.core = 3)
  DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
  DSM1 <- stack(paste("./DSM/",DSM[i],sep = ""))
  DSM0.Crop <- fieldCrop(mosaic = DSM0,fieldShape = Test.Crop,plot = F)
  DSM1.Crop <- fieldCrop(mosaic = DSM1,fieldShape = Test.Crop,plot = F)
  DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
  CHM <- DSM1.Crop-DSM0.R
  CHM.Rotate<-fieldRotate(CHM, theta = 2.3,plot = F)
  CHM.RemSoil <- fieldMask(CHM.Rotate, mask = EX1.RemSoil$mask,plot = F)
  EPH <- fieldInfo(CHM.RemSoil$newMosaic, 
                   fieldShape = EX1.Info$fieldShape, 
                   fun = "mean")
  DataTotal<-rbind(DataTotal,
                   data.frame(DAP=as.character(do.call(c,strsplit(MOSAIC[i],split = "_"))[2]),
                              EPH$fieldShape@data))
  # Making map plots
  fieldPlot(fieldShape=EPH$fieldShape,
            fieldAttribute="NGRDI", 
            mosaic=EX1.Rotate, 
            color=c("red","green"), 
            # min.lim = 0,
            # max.lim = 0.35,
            alpha = 0.5,
            round = 2)
  print(paste("### Completed: ", "Mosaic_",i," ###",sep=""))
  }

# Organizing the data table:
colnames(DataTotal)<-c(colnames(DataTotal)[-c(dim(DataTotal)[2])],"EPH") # layer=EPH
DataTotal<-DataTotal[,!colnames(DataTotal)%in%c("ID","ID.1","PlotName")] # Removing column 12 ("ID.1")
DataTotal

# Saving extracted data "DataTotal.csv":
write.csv(DataTotal,"DataTotal.csv",row.names = F,col.names = T)

• Graphics

DataTotal$Name<-as.factor(as.character(DataTotal$Name))
DataTotal$Row<-as.factor(as.character(DataTotal$Row))
DataTotal$Column<-as.factor(as.character(DataTotal$Column))
DataTotal$DAP<-as.numeric(as.character(DataTotal$DAP))
DataTotal$NGRDI<-as.numeric(as.character(DataTotal$NGRDI))
DataTotal$BGI<-as.numeric(as.character(DataTotal$BGI))
DataTotal$EPH<-as.numeric(as.character(DataTotal$EPH))

ggplot(DataTotal, aes(x = NGRDI,fill=as.factor(DAP))) +
  geom_density(alpha=.5,position = 'identity') +
  facet_wrap(~DAP,ncol = 1)+
  scale_fill_grey(start=1, end=0)+
  labs(y="#genotypes",x="NGRDI", fill="DAP") +
  theme_bw() 

• Linear Regression

DataTotal.Reg<-subset(DataTotal,DAP=="40")
DataTotal.Reg$Check<-as.character(DataTotal.Reg$Name)
DataTotal.Reg$Check[!DataTotal.Reg$Check%in%c("G43","G44","G45")]<-""

ggplot(DataTotal.Reg,aes(y=NGRDI, x=EPH)) + 
  geom_point() +
  geom_smooth(method=lm)+
  labs(y="EPH",x="NGRDI",fill="",alpha="")+
  geom_vline(aes(xintercept=0.02),col="red", linetype = 2, size=0.7) +
  theme_bw()+
  geom_text_repel(aes(label = Check),
                  size = 3.5, col="red",
                  box.padding   = 0.35, 
                  point.padding = 0.5,
                  segment.color = 'grey50')+
  theme(legend.position = "none",
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) 

• Heritability

DAP<-unique(DataTotal$DAP)

H2<-NULL
for(h in 1:length(DAP)){
  mod<-lmer(NGRDI~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
  H2.a<-c("NGRDI",DAP[h],as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
  
  mod<-lmer(EPH~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
  H2.b<-c("EPH",DAP[h],as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
  
  H2<-rbind.data.frame(H2,rbind(H2.a,H2.b))
}
colnames(H2)<-c("Trait","DAP","H2")
H2$H2<-as.numeric(as.character(H2$H2))
H2$DAP<-as.numeric(as.character(H2$DAP))

ggplot(H2,aes(x=as.factor(DAP),y=H2,fill=as.factor(DAP)))+
  geom_bar(stat="identity")+
  facet_wrap(~Trait)+
  scale_fill_grey(start=0.8, end=0.2)+
  labs(x="Days After Planting (DAP)", fill="")+
  geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
  theme_bw()+
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Area under the curve (AUC)

DataTotal1<-DataTotal[as.character(DataTotal$Name)%in%c("G43","G44","G45"),]

ggplot(data=DataTotal1, aes(x=as.numeric(DAP), y= NGRDI, col= Name, group=Name)) +
  geom_point(size=6)+
  geom_line(size=1.2) +
  scale_color_grey(start=0.8, end=0.2)+
  labs(x="Days After Planting (DAP)", fill="", col="")+
  theme_linedraw()

Trait<-c("NGRDI") # c("GLI","EPH")
Plot<-as.character(unique(DataTotal$Plot))

DataAUC<-NULL
for(a1 in 1:length(Plot)){
  D1<-DataTotal[as.character(DataTotal$Plot)==Plot[a1],]
  x1<-c(0,as.numeric(D1$DAP))
  y1<-c(0,as.numeric(D1[,Trait]))
  DataAUC <- rbind(DataAUC,
                   c(NGRDI_AUC=AUC(x = x1[!is.na(y1)], y = y1[!is.na(y1)]),
                     Name=unique(as.character(D1$Name)),
                     Trait=unique(D1$Trait),
                     Row=unique(D1$Row),
                     Column=unique(D1$Column)))}

DataAUC<-as.data.frame(DataAUC)
DataAUC$NGRDI_AUC<-as.numeric(as.character(DataAUC$NGRDI_AUC))
DataAUC$Name<-as.factor(DataAUC$Name)
DataAUC$Row<-as.factor(DataAUC$Row)
DataAUC$Column<-as.factor(DataAUC$Column)
DataAUC                    

• AUC Heritability

# Mixed model:
mod<-lmer(NGRDI_AUC~Row+Column+(1|Name),DataAUC)
H2<-as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov)
H2

ggplot(DataAUC, aes(x = NGRDI_AUC)) +
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  geom_density(alpha=.5,position = 'identity', fill="cadetblue") +
  labs(y="#genotypes",x=paste("Area under the curve (AUC_NGRDI: H2=",round(H2,2),")",sep="")) +
  theme_bw()